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We use Molecular Dynamics simulations to study attractive interactions and the underlying ionic 
correlations between parallel like-charged rods in the absence of additional salt. For a generic 
bulk system of rods we identify a reduction of short range repulsions as the origin of a negative 
osmotic coefficient. The counterions show signs of a weak three-dimensional order in the attractive 
regime only once the rod-imposed charge-inhomogeneities are divided out. We also treat the case 
of attraction between a single pair of rods for a few selected line charge densities and rod radii. 
Measurements of the individual contributions to the force between close rods are studied as a function 
of Bjerrum length. We find that even though the total force is always attractive at sufficiently high 
Bjcrrum length, the electrostatic contribution can ultimately become repulsive. We also measure 
azimuthal and longitudinal correlation functions to answer the question how condensed ions are 
distributed with respect to each other and to the neighboring rod. For instance, we show that the 
prevalent image of mutually interlocked ions is qualitatively correct, even though modifications due 
to thermal fluctuations are usually strong. 



I. INTRODUCTION 

The interaction of polyelectrolytes or charged colloids 
in polar solvent depends sensitively on the structure of 
the electrical double layer surrounding the macroion. 
This double layer consists in the simplest case of small 
counterions of opposite charge. The linearized mean-field 
treatment of this layer lies at the heart of DLVO theory 
[0, 1^ , one of the most influential and still very important 
descriptions of these systems. 

Today it is well established that sufficiently strong elec- 
trostatic interactions entail phenomena qualitatively be- 
yond the mean- field level. The possibility of attractive 
interactions between like-charged cylindrical macroions 
due to correlated ion fluctuation has been noted as early 
as 1968 by Oosawa [|[ ||, and Patey showed in 1980 
(using an integral equation approach and the HNC clo- 
sure) that two like-charged spheres will ultimately at- 
tract HJ. For various reasons these results were initially 
viewed with some scepticism, but attractive interactions 
and other non-mean-field phenomena (like overcharging) 
have soon after been confirmed by a large number of stud- 
ies, based for instance on computer simulations ^ ^, |^, 
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density functional theories 
theoretical calculations |^ , [13| 
proaches @ ||, || ||, If |2|, |5^, 
excellent recent summary can be found in Refs. |59 
Furthermore, experiments have shown that DNA (a stiff, 
highly negatively charged polyelectrolyte) can be con- 
densed by multivalent counterions 





This correlation-induced attraction is for instance be- 
lieved to be important for the compaction of DNA inside 
viral capsids |6^, ^ . 

Even though many of the abovementioned theories are 



quantitatively very successful, they are often also math- 
ematically fairly involved — like any systematic attempt 
to improve upon mean-field theory. It is thus desir- 
able to have a theoretical description which - based on 
some knowledge about the nature of the dominant cor- 
relations - gives a more direct insight into the physics. 
For example, the Wigner-crystal theories, starting with 
Refs. [^^l], ^ follow these lines and estimate the 
excess correlational free energy by looking at the ordered 
ground state of a two-dimensional layer of adsorbed ions. 
The success of these semi-empirical approaches clearly 
depends on how well one understands the existing corre- 
lations. However, even though the effects of correlations 
on, say, effective potentials or the phase behavior have 
been well studied in the past, their nature has attracted 
much less attention ^ |9). 

Recently several theoretical models have been pub- 
lished which discretize the distribution of condensed ions 
on DNA (by assuming occupiable lattice-sites) and treat 
the resulting partition function analytically or numeri- 
cally ^ |5^, Q . These works provided a further im- 
portant step in understanding the origin of attraction and 
the nature of the correlations involved, but it is not al- 
ways obvious how the employed discretization influences 
their strength. 

In the present paper we use Molecular Dynamics (MD) 
simulations in order to study the nature of counterion 
correlations around charged cylindrical polyelectrolytes 
(modeled as charged rods) in the absence of additional 
salt and on the level of a dielectric continuum approxi- 
mation for the solvent. We will investigate both the case 
of a bulk system of parallel rods, in which the occurrence 
of a negative pressure is a sure indicator of attractions, 
as well as a single pair of rods, in which the total force 
between the macroions is the appropriate observable. Re- 
call that we invariably are concerned with effective free 
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energies of interaction after "integrating out" ionic de- 
grees of freedom, which generaUy renders these interac- 
tions non pairwise additive. Hence, our two sets of data 
complement each other. In both cases it provides addi- 
tional insight to look separately at contributions coming 
from electrostatic and non-electrostatic origin. 

Finally we would like to remind the reader that short- 
ranged correlation induced attractions have to be care- 
fully distinguished from the much longer ranged attrac- 
tion between like-charged macroions, indications of which 
have been found experimentally in dilute suspensions of 
highly charged colloids at low ionic strength (for recent 
work see for instance Ref. g ^ ^, 0). Suggestions 
for a theoretical explanations of this phenomenon rest 
on reconsiderations of DLVO theory (7|, |7|, |7|, |7|l , but 
the situation is much less clear-cut — both experimentally 
IttI as well as theoretically In the present 

article we will not be concerned with these phenomena. 



II. THE SIMULATION METHOD 

Simulation details for the bulk systems have been de- 
scribed in Refs. |8|, In brief: We perform MD 
simulations using a Langevin thermostat to drive the sys- 
tem into the canonical state jsj]. We will need only two 
kinds of basic interactions. First, a short range repulsion 
which we model by the repulsive part of a Lennard- Jones 
(LJ) potential 
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r < 2i/V 
otherwise 



(1) 



where a is essentially the distance below which strong 
repulsion sets in; we will use it as our unit of length. The 
energy scale is set by e, but for a purely repulsive LJ- 
potential its precise value will not matter and we set it 
equal to the thermal energy kBT. Eqn. (|l|) is sometimes 
also referred to as the "Weeks- Chandler- Andersen po- 
tential" Second, the bare Coulomb potential Vc{r) 
between two charges zie and Z2e can be written as 



PeVcir) = ziZ2- 



(2) 



where /3 = l/k-gT and the Bjerrum length = 
/3e^/47reoer measures the coupling strength by specify- 
ing the distance at which two unit charges have interac- 
tion energy k^T . For instance, using the relative dielec- 
tric constant of water e-c ~ 80 and ambient temperature 
T « 300 K we have Ib ~ 7 A. Under periodic bound- 
ary conditions the total Coulomb energy is obtained by 
a sum over all pairs - including the images -, for which 
we use efhcient particle-mesh routines |8^ . 

Our system consists of immobile charged rods and 
mobile oppositely charged counterions. The rods are 
assembled from a string of negatively charged spheres 
( "monomers" ) of diameter a which sit on a line at a sepa- 
ration of 6 = 1.042 (7. We place one rod along the main di- 
agonal of the central simulation box which under periodic 



boundary conditions yields a hexagonal array of infinitely 
long rods. The counterions are also modeled as spheres 
of radius a. If there are N counterions of valence v, 
global charge neutrality requires ^/SL/b = Nv and thus 
an average counterion density of n = N/L^ = ^/S/vbL^. 
No simulations presented in this paper contain additional 
salt. Let us finally introduce the dimensionless charge 
parameter ^ = /b, which measures the number of unit 
charges along one rod per Bjerrum length. We briefly re- 
mind the reader of the concept of Manning condensation, 
stating that if > 1, a fraction 1 — of counterions 
will associate with the rod |Q. A precise meaning of 
"association" is provided by the Poisson-Boltzmann so- 
lution of the cylindrical cell model [ p9| pO| and has been 
carefully discussed in the literature ]21|, |91|, |9^ , |9^ , 

In a second set of simulations we place two rods par- 
allel to the 2:-axis of the cubic box. Their separation is 
small compared to the box length in order to decouple the 
pair of rods from their periodic images. In this case we 
used both three-dimensional periodicity as well as only 
one-dimensional periodicity along the directions of the 
rods. The latter case was handled by a one-dimensional 
version of an algorithm that has been termed "MMM" 
iQ . This is an alternative method for evaluating electro- 
static forces in three-dimensions based on a convergence 
factor approach of order 0{N\ogN), and which can be 
adapted straightforwardly to two- and one-dimensional 
periodic systems. We evaluate the forces pairwise using 
two different formulas, one very efficient for particle pairs 
with large distances in the non-periodic plane, the second 
formula, which is slower, for particle pairs with small dis- 
tance. The formulas are derived similarly to the formu- 
las presented in [p5, 96 1. The formula for distant particle 



pairs uses a Fourier transform, while the formula for near 
particles uses an expansion in polygamma functions. In 
this complementary approach we also used a continuous 
line charge on the rod axis instead of an aligned assembly 
of point charges. We found identical results when using 
the three-dimensional and the MMM-based approach and 
will thus in the following not distinguish between results 
originating from either one. 

We fixed the surface-to-surface distance of the two rods 
to 2 cr, which is a typical separation at which attractions 
can be expected at sufficiently high Bjerrum length. For 
these systems we considered three different rod diameters 
To, namely ro/a = 0.5, 2.0 and 6.271, which were imple- 
mented by the repulsive part of a cylindrical Lennard- 
Jones potential, which was shifted towards larger radii 
by replacing r in Eqn. (|^) by r — (ro — cr/2). The sys- 
tem 3 with the largest radius has a ratio of rod radius to 
charge spacing of 6.0, which is similar to that of DNA, 
where it is 10 A/ 1.7 A « 5.9, hence we refer to it as the 
DNA-like system. The system parameters can be found 
in Table | 

Let us finally remark that the assumption of a homo- 
geneously charged rod or linearly aligned point charges 
is often only a first approximation, since many impor- 
tant stiff polyelectrolytes {e.g., DNA or actin filaments) 
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system 


ro/cr 


b/a 


s/a 


L/a 


1 


0.50 


1.042 


2.95 


134.5 


2 


2.00 


0.1538 


6.00 


42.00 


3 


6.27 


1.042 


14.54 


134.5 



TABLE I: Geometry of the simulated two-rod-systems, ro is 
the rod radius, b the separation of unit charges along the rod, 
s the axial separation of the rods and L the box length. The 
surface-to-surface distance of the rods is 1.95o" in case 1 and 
2 a in the other two cases. The ions are always trivalent. 



feature a helical charge distribution. Implications of this 
additional structure on many aspects of helix-helix inter- 
actions are discussed in detail in a series of theoretical 
papers by Kornyshev and Leikin ^ (see also the 
simulations in Ref. |2^). Since in the present work 
we are concerned with more generic questions about the 
origin of correlations, we will neglect this complication. 



III. BULK SYSTEMS 
A. Osmotic coefficient 

The osmotic coefficient p is defined as the ratio be- 
tween the actual pressure and a fictitious pressure that 
would act if all interactions were switched off ("ideal 
gas"). For instance, the osmotic pressure of a dilute so- 
lution of charged rods is overwhelmingly dominated by 
the counterions, but a certain fraction of them may be 
sufficiently strongly localized by the macroions such that 
they do not contribute their full share (one says they 
are "osmotically inactive"). In the regime of counterion 
condensation, > 1, the osmotic coefficient should ap- 
proach l/2^u < 1 in the infinite dilution limit, while for 
finite concentrations p is larger |9^ . 

Figure |l] shows our MD-results for the osmotic coeffi- 
cient of bulk systems which are characterized hy / a — 
1, h/a = 1.042, ^b/c = 1, and w = 3. The pressure is 
identified with the component of the stress tensor per- 
pendicular to the rods ||82| , and its contributions coming 
from electrostatic and non-electrostatic {i.e., entropic 
and excluded volume) origin are plotted separately. It 
can be seen that within a density range from roughly 
n = 8 X 10~'^cr~'^ to 6.5 x lQ~^a~^ , corresponding to rod 
separations between 6.8 a and 2.4 cr, the osmotic pres- 
sure is negative. If the rods were not forced to remain at 
fixed separations, they would phase separate, producing 
a condensate and a dilute phase. 

Surprisingly, the electrostatic contribution (which is 
always negative and always favors contraction) shows no 
particular features in the regime where p < Q. It rather 
appears that the negative pressure originates from a drop 
in the repulsive excluded volume contribution which is 
large enough such that the electrostatic contribution can 
win against the short range repulsion. This finding is 
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FIG. 1: Contributions to the osmotic coefficient p as a func- 
tion of counterion density n for a bulk systems of paral- 
lel charged rods characterized by ro/a = 1, b/a = 1.042, 
is/a = 1, and u = 3. The heavy dots and crosses repre- 
sent the non-electrostatic contribution (kinetic and excluded 
volume) and the negative of the electrostatic contribution, re- 
spectively (the curves are guides to the eye). The inset shows 
the total osmotic coefficient, and the limiting value of infinite 
dilution 1/2^d « 0.174 is indicated by an arrow. 



consistent with earlier simulation results |10 which used 
finitely replicated cells. A possible explanation of this 
phenomenon rests on the following tempting picture: The 
excluded volume contribution to the pressure stems from 
the force that ions exert on the oppositely charged rods. 
If the density increases, the rods come closer to each 
other and start to pull the ions away from their neigh- 
bors, thereby reducing this force. However, at too high 
concentrations rods will again repel by pushing onto each 
other via the ions in between. We will come back to this 



effect in Sec. IV B 



An alternative explanation suggests that the con- 
densed ions between two neighboring rods form a mutu- 
ally interlocked pattern, which results in a lower Coulomb 
energy (compared to a homogeneous charge distribution) 
and which thus leads to attractions ||l^, ^ |5^, |9| . Since 
each rod has six neighbors, this defines six "planes" in 
which one would have to look for a two-dimensional in- 
terlocking pattern. However, we were unable to find a 
corresponding structure in our simulated data that goes 
beyond a weakly developed first correlation hole, which 
on its own is not a sure sign for attractions as we will 
show in Sec. IV C for the case of two rods. This may 



partly be related to the small degree of arbitrariness in 
the definition of the planes, but there is also a physical 
explanation to consider. The systems are comparatively 
dense in the regime where attractions occur. It is there- 
fore likely that close ions are correlated irrespective of the 
plane they "belong" to — in other words, that those arti- 
ficial planes are irrelevant for understanding the physics. 
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It rather seems more reasonable to study the full three- 
dimensional correlations of the ions, as we will do in the 
following section. 



B. Pair correlation function 

If the radius of the rods and their mutual distance is 
smaller than the Bjerrum length, the electrostatic inter- 
action of ions reaches beyond the nearest rod and the 
complete system of ions may form a three-dimensional 
correlated structure. More specifically, Shklovskii sug- 
gests 1^ that in the case ro/h <^ v, or equivalently 
fo/^B <C w/C, ions could form such a structure on the 
background of the rods, which for very high coupling can 
be described as a three-dimensional Wigner crystal (or 
at least as a strongly correlated hquid). We tested this 
assumption by computing the usual pair correlation func- 
tion g{r) for the systems above which showed negative 
pressure. Before we discuss the results, a few general 
remarks are appropriate. 

If the ions form a three-dimensional crystal structure, 
there will be a typical minimum distance between them 
which should scale like one over the third root of the 



counterion density 111 ]. However, there is a second im- 
portant length scale, which is the separation of the rods 
and which is easily seen to scale as one over the square 
root of the counterion density. Hence, whatever the de- 
tails of the actual ionic structure are, if one measures a 
characteristic length varying like it can be viewed 

as a bulk correlated liquid property, while a length vary- 
ing like n~^/^ is rod- imposed. 

In Fig. H we present results for the system studied 
above in the density range 0.02 cr"^ ... 0.1 cr"^, which 
corresponds to the high-density range in which the os- 
motic pressure was found to be negative. The dashed 
curve shows the pair correlation g(r) at the density 
n = 4.25 X 10~^(T^^, which is at the minimum of the 
osmotic coefhcient in Fig. |l|. A pronounced oscillatory 
structure is clearly visible. The inset shows the position 
of the first maximum as a function of density in a double 
logarithmic plot (open circles). In this range the data 
suggest a power law with exponent —0.48, which is close 
to the value expected for a rod-imposed structure. In 
fact, the oscillations in g{r) simply reflect the periodicity 
of the array of rods. 

In order to extract actual inter-ionic correlations from 
the simulations, the imposed periodic inhomogeneity of 
the ion density has to be removed. One way of doing this 
is as follows: The pair correlation function is the proba- 
bility of finding an ion at a distance r from another ion 
relative to the probability of finding an ion at the same 
distance in a random (i.e., noninteracting, "ideal gas") 
system at the same average density. In the present case 
it would be sensible to normalize not by a random homo- 
geneous system but by a random inhomogeneous system, 
i.e., a system in which there are no inter-ionic correla- 
tions but the spatially varying ion density is preserved. 
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FIG. 2: Pair correlation function g(r) for the bulk system 
from Fig. |l| at a counterion density n = 4.25 x 10~^cr~^ 
(dashed line). The solid line is the ratio between this g(r) 
and grand (f) (see text). The inset shows a double logarithmic 
plot of the position of the first maximum of g(r) (open cir- 
cles) and (?(f')/grand('") (fuU circles) as a function of density; 
the slopes are indicated. 



This can be easily accomplished in the following way: In 
each configuration move every ion a random distance (be- 
tween and \f2)V) along the direction of the rods. This 
does not change the inhomogeneous one-particle distri- 
bution, but completely destroys all two-particle corre- 
lations. Then compute the pair correlation function of 
this randomized system, grand (?")• The ratio between the 
usual g(r) and the randomized (7rand('') now contains all 
information about inter-ionic correlations, but the im- 
posed density inhomogeneity is divided out. This ratio 
is also plotted in Fig. ^ for the system at the density 
n = 4.25 X 10~^(T~^ (solid line). Signs of correlations can 
again be seen, but the long range oscillatory part has 
been removed. The inset also shows the position of the 
first maximum of gij) j g■ca.r^A{r) as a function of density 
(closed circles). The solid line has a slope of —0.35, which 
is more close to the exponent expected for a correlated 
three-dimensional ionic structure. 

We see that the nature of correlations in this sys- 
tem is a subtle interplay between rod-imposed and inter- 
ionic contributions, and the latter are only identifiable 
once the former are divided out. The idea that a three- 
dimensional pattern of correlated ions forms on the struc- 
tureless background of the rods is clearly too simple for 
the system under study. Furthermore, both correlations 
are comparatively weak, in the sense that the first max- 
imum in g(r) is fairly low, and the long range structure 
visible in the bare g(r) is identical in grand (?") and thus 
is imposed by the external periodicity of the rods. The 
observed correlations are definitely much weaker than ex- 
pected for a Wigner crystal or a strongly correlated liq- 
uid. If the ideas advocated in Ref. ||5ll, |5^, ^ about 
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the source of attraction are correct, the present analysis 
impUes that a very low degree of correlations is already 
sufhcient or, in other words, the Wigner crystal picture 
holds way beyond the ground state. This of course pro- 
vokes the questions "Why?" and "How far beyond?" . 

The above analysis does not directly explain our earlier 
finding that the occurrence of a negative pressure is ulti- 
mately related to a sudden drop in the repulsive excluded 
volume forces. One may speculate that the electrostatic 
part of the pressure responds less sensitively than the 
short range LJ-part to a comparatively local ordering 
as observed in Fig. ^, but additional studies would be 
needed to support this. 



IV. ONE PAIR OF RODS 

The Wigner crystal picture is ultimately based on en- 
ergetic arguments, it does not provide a direct "mechan- 
ical" explanation for why correlations actually produce 
an attractive force. To answer this question, these corre- 
lations have to be studied in more detail. However, the 
bulk system is not necessarily the easiest case: The rele- 
vant observable is the pressure, which is less direct than 
a force, and the more complex geometry complicates the 
definition of suitable and reasonably intuitive correlation 
functions. We therefore resort now to the interaction be- 
tween two rods as well as the involved ionic correlations. 
As we have mentioned in the introduction, studying the 
pair interactions is not merely an alternative but also a 
complementary approach to studying the bulk system, 
since the forces in the latter cannot be decomposed into 
pair forces. 

A pair of parallel rods has previously been investigated 
as a function of rod separation, which gives the distance 
dependence of the force |l3 , 20 , 5^, |5^ and by integration 
the potential of mean force. In our work we supplement 
these results by studying the properties of the system at 
fixed distance as a function of Bjerrum length. This al- 
ternative scan is promising since unlike for bare charges 
the Bjerrum length does not simply enter as a prefac- 
tor of the interaction. Rather, its size will determine to 
what degree ions can be found close to the rods and how 
strongly they correlate with it as well as with each other, 
thereby infiuencing the nature of the force itself. 

We studied three systems which differ in their values 
for rod radius and line charge density, as specified in 
Tab. |, but we always kept the surface-to-surface dis- 
tance between two rods around 2 a, i. e., twice the diam- 
eter of a counterions, since this is a typical separation 
at which attractive forces have been found at sufficiently 
large Bjerrum length. The observables we measured were 
the force between the rods (split again into electrostatic 
and non-electrostatic origin) as well as various ionic cor- 
relation functions along and around the rods. The three 
studied systems will provide examples for qualitatively 
different behavior, but cannot predict detailed dependen- 
cies on key system parameters. For instance, the simula- 



tions in Ref. suggest that increasing the rod radius at 
fixed line charge density will ultimately remove repulsion, 
while increasing the line charge density at fixed rod ra- 
dius will ultimately lead to attraction. However, present 
stud ies in this direction show a more complex behavior 

It is well known that the importance of correlations 
is strongly influenced by the valence of the counterions. 
However, this variable cannot be changed continuously. 
In contrast, increasing the Bjerrum length can be viewed 
as continuously increasing the strength of correlations 
and thereby studying in detail how they develop. Experi- 
mentally the Bjerrum length is a parameter whose change 
requires some effort, but it can be achieved for instance 
by a careful control of the solvent dielectric constant, as 
has been recently shown in an experimental study of the 
coil-globule transition of DNA |101|. 



A. Force as a function of Bjerrum length 

In a first step we study the force between the two rods 
by (separately) adding up all electrostatic and excluded 
volume forces that act on one rod and average over the 
simulation run (which typically means over about 3000 
configurations separated by 2000 integration steps). Un- 
der periodic boundary conditions the total force on one 
rod originates not just from the neighboring rod, but also 
from all periodic images. However, we assume these con- 
tributions to be negligible in our case for the following 
reasons: For high Bjerrum length all ions are condensed 
onto the surface of the rods, implying that interactions 
of these essentially neutral objects {i. e., rods plus coun- 
terions) over a distance of the box length is much weaker 
than the direct interaction (compare the values for s and 
L in Tab. |[). It is not immediately clear that this re- 
mains valid in the limit ^ 0, but in this case our data 
are always asymptotic to the analytical expression for the 
electrostatic force between two charged lines (see below) , 
which implies that the images can be neglected as well. 

One may ask why we used periodic boundary condi- 
tions in the first place, if we want images to be irrelevant. 
The answer is of course that the images along the direc- 
tion of the rod are important, since they render the rods 
infinitely long — or stated differently: remove end effects. 
Periodic replication perpendicular to the axis of the rods 
is an unwanted but acceptable (since controllable) side 
effect imposed by the use of traditional Ewald schemes, 
in which it is not trivial to switch off periodicity in a 
specific direction. However, the latter can be achieved 
quite easily by using a one-dimensional version of MMM 
p\ |9^ , p6| , compare Sec. II, and we will also present data 
obtained by this second method. 

Motivated by our observations in bulk simulations we 
also look at the components of the force originating from 
electrostatic and excluded volume interactions. Let us 
start with a few general remarks concerning what to ex- 
pect. For €b = there is no electrostatic interaction 
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FIG. 3: Force between two rods of system 1 from Tab. | as a 
function of Bjerrum length £b- Total force, excluded volume 
contribution and electrostatic contribution are represented by 
"•" on a solid line, on a dashed line, and "x" on a dotted 
line, respectively. The fine dashed line indicates the limit- 
ing behavior from Eqn. (^. Positive forces denote repulsion, 
negative attraction. Measured values are indicated by the 
symbols, the lines should merely guide the eye. 



between the rods, and the onl y pos sible source of inter- 
action is a depletion attraction |102] generated by the ex- 
clude d vo lume of the (at = effectively uncharged!) 
ions 1 112]. When slightly increasing the Bjerrum length, 
the two rods will feel an unscreened Coulomb repulsion, 
since for < 1/2 ions will remain unbound |p8| . The 
force per unit length is then given by 



0F ^^=" 



(3) 



i. e., proportional to Bjerrum length. In other words, in 
this regime the Bjerrum length still acts as a prefactor de- 
termining the strength of the interaction. Upon further 
increase of ion condensation and subsequent correla- 
tions set in, resulting in nontrivial i^-^B-curves (see be- 
low). However, in the limit £b — oo everything becomes 
simple again: The counterions will ultimately assume a 
"ground state" configuration which will no longer change 
upon further increase of £b (only fluctuations around the 
ground state will become weaker). Hence, the electro- 
static force (and as an indirect consequence also the ex- 
cluded volume force) will again be proportional to £b- 
However, the constant of proportionality cannot be pre- 
dicted from these considerations — not even its sign. In 
the following we will plot the force per unit length and 
use the convention that a positive sign denotes repulsion^ 
while a negative sign denotes attraction. 

Fig. H shows the result for system 1 (for notation see 
Tab. |), which consists of rods having the same diameter 
as the ions and the same separation as in the system from 
Fig. showing the most negative pressure. For small 
Bjerrum length the force is given by Eqn. (0). However, 



if > 1/2 (i.e., ^b > 0.17(t the present case) Manning 
condensation will set in on the two-rod system. Since ions 
will condense preferentially between the rods, they reduce 
the electrostatic repulsion, but at the same time produce 
an outward pressure due to their excluded volume. 

Provided the excluded volume contribution to the in- 
teraction is not yet substantial, a description of the sys- 
tem using Poisson-Boltzmann (PB) theory can work up 
to this Bjerrum length. Indeed, the PB equation can be 
solved exactly in this geometry ]103|, 104 , but the analyt- 



ical solution (in terms of hypergeometric functions) is ex- 
tremely complicated. One particularly direct result how- 
ever (exact for infinitely thin rods) is that if the two-rod- 
system is below the Manning threshold, pure Coulomb 
repulsion [i.e., unmodified by the presence of counteri- 
ons) is found asymptotically at large separation, while if 
each rod is at the Manning threshold, the Coulomb re- 
pulsion at large distances is reduced by a factor of 2. We 
mention aside that for the case of two rods (of arbitrary 
radius) and added salt the PB equation has been solved 
numerically [105|, and analytical studies i n th e presence 
of salt exist within linearized PB theory |106], even for 
tilted rods pO^ . 

Upon further increase in Bjerrum length the electro- 
static force is seen to change sign at ^b/o" ~ 0.6, and the 
total force becomes attractive beyond Ib/ct ~ 0.9. While 
it is easy to imagine counterion distributions between 
the two rods that would lead to electrostatic attraction, 
those distributions are counteracted by both entropy and 
excluded volume interactions among the ions, since a 
high density between the rods is req uired . A mean -field 
treatment on the level of PB theory lOj , 104 , 105 can- 
not resolve the issue because rigorous proves ex i st th at 



the above situation must give repulsion |108, 109, 110 



contrary to the actual observation, which is in agreement 
with earlier simulational works (l2[ |l^. It is clear that 
(on the level of the restricted primitive model of elec- 
trolytes) the attraction must hence be related to the ex- 
istence of correlations between the ions. We will discuss 
a few of them in the following section. 

For a Bjerrum length larger than 2.5 a the excluded 
volume part of the force becomes attractive as well. This 
is surprising, since it implies that there are still suffi- 
ciently many ions between the rods to induce electro- 
static attraction, but they organize their positions such 
that they put less pressure on the rods than the ions 
located on the outward surfaces. 

In Fig. ^ we show the same Bjerrum-scan for system 
2, which differs from the previous one in the following 
ways: The linear charge density is 6.78 times larger, the 
rod radius is four times larger, and the rods are kept at 
a separation of 6 a, which ensures that the surface-to- 
surface separation is again 2 a. The generic behavior at 
small Bjerrum length is the same, but in this system a 
very pronounced difference occurs at larger values: Be- 
yond £b/<^ ~ 1-5 the electrostatic contribution to the 
force becomes repulsive and the total force is attractive 
only because of the excluded volume term. Note that to- 
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FIG. 4: Force between two rods of system 2 as a function of 
Bjerrum length. The inset magnifies the initial region Is £ 
[0; 0.6]. The line styles are the same as in Fig. ^ 
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FIG. 6: Relative imbalance SNin = (Mn - A^out)/(A''m + A^'out) 
of ions between the rods and outside (see text) as a function 
of Bjerrum length for the three systems studied. Note that a 
positive SNin means that more ions are between the rods. 




FIG. 5: Force between two rods of system 3 as a function of 
Bjerrum length. The inset magnifies the initial region £b £ 
[0; 3]. The line styles are the same as in Fig. ^. 



gether with the results discussed above this implies that 
a net attraction can occur both because electrostatic at- 
traction overcomes excluded volume repulsion and be- 
cause excluded volume interactions overcome an electro- 
static repulsion. In the present case electrostatic repul- 
sion is only weakly developed and numerical errors are 
large compared to the other two systems, but we have 
observed the same phenomenon in several other systems 
as well. Its characteris tics will be discussed more thor- 
oughly somewhere else [ |100| . Even though this effect may 
appear counterintuitive, we want to remind the reader 
that minimizing the total energy of the ionic system at 
given rod separation does not imply that the electrostatic 
part of it gives rise to rod-rod-attractions. 



Fig. H presents results for system 3, which compared 
to system 1 has an approximately 12.5 times larger rod 
radius. The surface-to-surface separation is again main- 
tained at 2 a. The ratio ro/6 « 6 is close to that for DNA. 
Fixing the length scale via b = &dna = 1.7 A implies a 
(somewhat small) ion diameter of cr = 1.63A Im|. At 
the Bjerrum length = 7.14 A = 4.38 cr (as appropriate 
for water at room temperature) the total force is attrac- 
tive, and indeed it has been experimentally established 
that trivalent ions lead to attractive interactions between 
DNA strands [|l], |62[ H |6|, |65[ H . Note also that just 
like in system 1 both contributions to the force are at- 
tractive at sufficiently high Bjerrum length; however, in 
system 3 the excluded volume part is stronger than the 
electrostatic part for values of ^b/c larger than « 6. 

As we have mentioned above, the attraction between 
the rods due to electrostatic forces arises because ions 
condense preferentially between the rods. A straight- 
forward measure for this imbalance is provided by the 
following observable: Imagine two parallel planes, each 
containing the axis of one of the rods, whose distance 
equals the axial distance between the rods. These planes 
divide space into a region "between" the rods and two 
disjunct regions "outside" . Denote by iVi„ and A^out the 
number of counterions in the region "between" and "out- 
side" , respectively, and define the relative ionic excess 
between the rods as SN-^ = (iVi„ - Nont)/{Nin + Nout)- 
This observable is plotted in Fig. ^ as a function of Bjer- 
rum length for systems 1, 2 and 3. For small Bjerrum 
length SN\a is negative, since ions are not condensed and 
the outside region is larger. As ion condensation sets in 
at increasing £b, SNin also increases until it reaches 0, 
the point at which the inside and outside numbers are 
balanced. For system 1, SN^n becomes strongly positive 
afterwards, showing that there are more ions between 
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FIG. 7: Geometry for the definition of the azimuthal corre- 
lation function g[Lp\5r). It is proportional to the probability 
density of finding a counterion at the angle <^ relative to the 
other rod within a condensation distance of at most Sr. The 
normalization is such that g = 1 corresponds to the average 
density, or stated differently, the average of g is 1. 



FIG. 8: Polar plot of azimuthal correlation functions g{(p\s/2) 
for system 1. The particular Bjerrum lengths are indicated. 
The neighboring rod is assumed to be located to the right 
side. The normalization is such that the average value of g is 
equal to 1 (indicated by the dotted circle). Measured values 
(dots) are only indicated for the case Ib ~ lo". 



the rods than outside, but it decreases again beyond 
i^/a ~ 1.5, implying that this imbalance softens out. 
For system 2, SN-m does not rise significantly above 0, but 
rather approaches this balance line at about £b/o' ~ 1. 
For system 3, 6Nin weakly rises above at ^b/ct ~ 1.6, 
but it again drops below at ^b/c ~ 3.6. 

In all three systems the point at which the electro- 
static force becomes attractive coincides roughly with the 
point at which a strongly negative imbalance SN-^ van- 
ishes. However, one has to be careful in interpreting this 
observable: For system 1 the excluded volume repulsion 
between the rods vanishes beyond ^b/c ~ 2.5, but 6Nin 
is positive there, i.e., there are still more ions between 
the rods than outside. The attraction between the rods 
in system 2 is due to excluded volume forces, but the 
imbalance 6Nin is essentially 0. Hence, attraction occurs 
not due to a simple density difference. In the following 
section we will therefore also study the azimuthal distri- 
bution of ions in some more detail as well as the question 
how close the ions actually come to the rods. 

Let us close our discussion of forces with the following 
observation. Compared to system 1 the density of sur- 
face charges ^ is 12.5 times lower in system 3, but the 
total force becomes attractive at ^b/c ~ 2.9, which is 
only 3.2 times larger than for system 1. Rouzina and 
Bloomfield [ [49| use the two-dimensional plasma parame- 
ter r2 = fBi'^^T^^^ as a measure for the strength of cor- 
relations and estimate that the onset of attraction should 
be expected at r2 « 2. This implies that the correspond- 
ing Bjerrum length scales inversely proportional to the 
square root of the surface charge density. For our simu- 
lations r2 ~ 2.58, 1.83 and 2.35 for systems 1, 2 and 3 
at the onset of attraction, which agrees remarkably well 
with their estimate. 



B. Azimuthal correlations 

The previous sections have revisited the importance of 
ionic correlations for the nature of the force between two 
charged rods, and in the following we will present mea- 



surements of a few of them. We start with an ion-rod 
correlation that answers the question how the ions dis- 
tribute around a rod relative to the position of the other 
rod. Fig. illustrates our definition of the correlation 
function. Denote the position of a counterion relative to 
rod 1 in polar coordinates r and ip, such that ip = corre- 
sponds to the direction towards rod 2. We define g{(p\6r) 
to be 27r times the probability density of finding an ion 
at the angle (p, given that its separation from rod 1 is at 
most Sr. Note that this implies that the mean value of 
g{ip\dr) (when averaged over ip) is 1. Obviously, g{ip\Sr) 
is periodic in ip with period 27r. At first we will choose 
Sr as large as possible, i.e., half the rod separation s. 
We want to mention that for Sr —^ ro this function is 
essentially the rod-analog of the density resolved by the 
polar angle for the spherical case studied in Ref. p5| |. 

Figure || shows a polar plot of g{(p\s/2) for system 1. 
The Bjerrum length is varied between 1 a and 5 a. Com- 
pare also the corresponding force-plot from Fig. ||. Sev- 
eral things can be observed: (z) The azimuthal correla- 
tion function is largest in the direction towards the other 
rod. (a) There exists a "second tide" on the side oppos- 
ing the other rod. (Hi) The elongation of f/((p|s/2) along 
the line joining the two rods increases with increasing 
Bjerrum length, (iv) The larger than average g at ip = 
and maybe = tt is balanced by a less than average g 
in the transverse direction, (v) While g{0\s/2) is always 
larger than 1, g{TT\s/2) is smaller than 1 for sufficiently 
small Bjerrum length, but ultimately rises above 1 as 
well. 

As can be seen from Fig. ||, all five systems presented in 
Fig. H feature a net attraction, but for the systems with 
^b/c = 1 and 2 the excluded volume part is still repul- 
sive. It is unfortunately difficult to directly relate this 
finding to the shape of the azimuthal correlation func- 
tion, since the position of the ions are of course relevant 
for both electrostatic and excluded volume forces, which 
are opposite and show a different distance dependence. 
We come back to this issue below, when we discuss the 
dependence on dr. 
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FIG. 9: Same plot as in Fig. |8| for system 3. Measured values 
(dots) are only indicated for the case £b = Icr. 



Fig. ^ shows exactly the same five correlation functions 
for system 3. This system differs from system 1 "only" 
by a 12.5 times larger rod radius, but the shape of the 
correlation functions is quite different. The second peak 
at (/? = TT has given way to two new peaks emerging some- 
what beyond 7r/2, while the ion density at (/? = tt is below 
average. The fact that the peak at (/? = decreases upon 
increasing Bjerrum length can be traced back to the nor- 
malization of g and the fact that the total number of 
ions within the shell r < s/2 initially increases as £b in- 
creases. At low Bjerrum length most of the condensed 
ions can be found between the rods. Upon increasing 
£b the additionally condensed ions will occupy the other 
parts of the rod, thereby reducing the relative weight of 
the ones at = 0. Once essentially all ions have con- 
densed within the shell r < s/2, a larger Bjerrum length 
will only enhance the correlations and thereby make the 
peak at = larger again. We note that in the spirit of 
our mapping indicated above the azimuthal correlation 
function for ions around DNA would approximately be 
given by the curves for ^e/f = 4 or 5. 

In Fig. ^ we show g{ip\s/2) for two versions of sys- 
tem 2. This system has a rod radius 4 times larger 
than system 1 and a line charge density about 6.8 times 
larger. For a Bjerrum length of — 0.5 a, for which 
the total force between the rods approximately vanishes 
(see Fig. a pronounced peak in g((p\s/2) towards the 
neighboring rod as well as a slightly smaller than av- 
erage counterion density on the opposite side is visible. 
Just as for system 3, upon increasing the Bjerrum length 
g{ip\s/2) does not merely "intensify" its features but de- 
velops a qualitatively new structure. In addition to the 
main correlation peak at f — 0, five new peaks show 
up, roughly at multiples of 60°. This occurrence of ad- 
ditional peaks is a common phenomenon in usual liquids 
as correlations increase. However, in the present case we 
have the additional constraint that g{ip\5r) is periodic in 
if and obviously an even function, which restricts the lo- 
cations of these peaks. While for a flat two-dimensional 
system the position of peaks in the pair correlation func- 
tion is solely determined by the density of charges (<j~^/^ 
is the only available length scale), in a curved geometry 
two more effects play a role: first, the radius of curvature 
appears as a new length scale, and second, topological 
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FIG. 10: Same plot as in Fig. |8|for system 2. Measured values 
for 1^1 a = 0.5 and Is/o- = 3.0 are indicated as open circles 
(o) and dots (•), respectively. The lines should merely guide 
the eye. 
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FIG. 11: Azimuthal correlation functions g{tp\5r) for the sys- 
tem 2 with £b = 3a (see also Fig. |l^ for different values 
of Sr. Note in particular that for decreasing Sr the peak at 
if = turns into a correlation hole. 



constraints require the function to close upon itself as 
in the case above. Particularly the latter observation 
suggests that there are situations in which this match- 
ing works "automatically" whereas in other situations it 
leads to a frustration. This may favor counterion confor- 
mations which relax these constraints {e.g. helices) and 
influence the strength of the rod-rod-interaction. 

We conclude our discussion of the azimuthal correla- 
tion function by discussing its dependence on Sr. This 
variable determines which counterions are taken into ac- 
count for computing the density as a function of if. If Sr 
gets smaller, the focus is on ions that are more closely in 
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contact with the rod. Fig. ^ shows a plot of g{ip\5r) for 
system 2 with f b = 3 it. The bold solid line is the same 
curve as shown in Fig. 10 (i.e., Sr = s/2 ~ Su), only 
plotted in a Cartesian way. The other three curves cor- 
respond to successively smaller values of 6r. We would 
like to draw attention to the fact that the peak at ip = 
gives way to a correlation hole. While g{Q\Sr) decreases 
upon reduction of Sr, g{T:\Sr) increases. This is impor- 
tant since the ions most closely in contact with the rod 
produce the strongest excluded volume force. Fig. |ll| 
thus shows that ions very close to the surface of the rods 
are more likely to push them inwards rather than out- 
wards. Note that this particular system is special in 
that the attractive forces have been found to originate 
from the excluded volume term, see Fig. ^. A similar 
(5r-analysis for system 1 in the strongly attractive regime 
at £b — "20(7 does not show this effect (data not shown) . 
As 6r gets smaller than 1.0 cr, the peaks at ip = and 
if = TT start decreasing and rising, respectively, but the 
former does not turn into a correlation hole. The same 
finding applies to system 3 in the attractive regime at 
£b = 5 cr. However, note also that in the latter two sys- 
tems electrostatics is the dominant reason for attraction. 
Unfortunately, effects at small 6r are difhcult to observe, 
since the small number of ions at these close distances 
impedes a good statistics. 

We want to remark that a similar correlation hole has 
been previously observed in the spherical case and has 
been termed an electrostatic depletion effect |15 . Ob- 
serve that in our case this effect is only visible if one fo- 
cuses on ions very close to the surface of the macroion — 
the total amount of ions between the macroions is above 
average, which contradicts the depletion picture. 



C. Ion interlocking 

It has previously been suggested that the ionic cor- 
relations in the case of attracting rods take the form 
of an interlocked pattern in the plane between the rods 
p5| , p9|| . In this section we present first measure- 
ments of these kind of correlations under off-lattice con- 
ditions. It turns out that alternating charge-asymmetries 
indeed exist, but they are remarkably weakly pronounced 
and fairly short ranged. We also investigate their rele- 
vance for the attraction. 

Fig. envisages a tentative "ground state" for system 
1 from Tab. |. We have previously deduced from Fig. ^ 
that in the limit of high Bjerrum length the ions are es- 
sentially either at cp = or ip = n. On each rod, every 
three monomers a trivalent counterion is condensed. The 
ions between the rods form an interlocked pattern that 
also gets imprinted on the ions located on the opposite 
side of the rods. Such a pattern suggests the definition 
of the following three pair correlation functions, as indi- 
cated in Fig. |l^: Given an ion, that is condensed on one 
rod and sits between the two rods, what is the probabil- 
ity of having a second ion at a distance z along the rods 
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FIG. 12: Scheme of a possible "ground state" of the two rod 
systems, showing ion interlocking. The relative position of 
three kinds of condensed counterions with respect to a con- 
densed ion between the rods give rise to three different corre- 
lation functions, as discussed in the text. 



3.0 



-1 1 1 r— 



-1 1 1 1— 

1.3 



— 1 1 1 r- 




FIG. 13: Interlocking correlation functions gi{z) (solid line 
g2{z) (dashed line) and gsiz) (dotted line) as defined in Fig. h 
and the text for system 1 with Bjerrum length £b = 5 a. The 
inset shows the ratio g2(z) / g-i[z). 



that is condensed (1) on the opposing rod and also be- 
tween the rods, (2) on the same rod and also between the 
rods, and (3) on the same rod, but facing outward? With 
the usual normalization to 1 at large distance, let us call 
these three functions gi(z), g2(z) and 53(2), respectively. 
We only consider ions whose center has a distance of less 
than cr/2 from the common plane of the rods. 

Fig. |l3| shows the result of a measurement of these 
three correlation functions for system 1 with = 5 a. 
Pronounced correlations are clearly visible, but their 
range is comparatively short. Structural features in the 
correlation function 53 are much weaker developed. The 
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reason for this is that the ionic density outside the rods is 
lower than between them, and the structure arises from 
a weak "transfer" of the correlations of gi and especially 

52- 

It may seem surprising that not just gi{z) but also 
g2iz) has its first peak at z « 36. The reason for this 
is that just as for the three -dime nsional pair correlation 
function discussed in Sec. [II B a fair amount of long- 



rangcd periodicity in the gi{z) arises due to an external 
imprinting. Here it is the fact that charge neutrality 
essentially requires 1 trivalent ion per rod every 3&. A 
visual inspection of typical ion configurations also shows 
that the system is still far from a highly ordered ground 
state as envisaged in Fig. ^ In particular, for £b/<^ — 5 
the imbalance SN-m presented in Fig. ^ has a value of 
about 0.022, which implies that the average density be- 
tween the rods is about 5% higher than outside. These 
additional ions have to be put somewhere and clearly 
destroy the simple groundstate from Fig. |l^. Together 
these findings imply that 3& is also a typical distance 
that has to be expected for 52, but observe also that the 
peak in 52 (-s) at z ~ 3& is less pronounced than the one 
for gi{z). In fact, the terminology "interlocking" sug- 
gests that we should be interested in the relative charge 
asymmetries along the rods. A better measure would 
therefore be the ratio between the two inner correlation 
functions, as is shown in the inset of Fig. For z ^ 
the ratio g2{z) / gi{z) approaches 0, since the correlation 
hole of g2{z) at z = is of course more pronounced than 
that for gi{z). The small peak shortly before z = 3b in- 
dicates the high probability for this distance as discussed 
above, but the fact that the peak value is below 1 means 
that such a value occurs more likely for ions belonging 
to different rods. A somewhat broader peak can be seen 
around z « 6b. It has its value above 1, indicating that 
this distance occurs preferentially for ions on the same 
rod. Beyond z = lOb no significant structure is visible. 

Together these findings show that ion interlocking does 
indeed occur, but that it is again superimposed by "triv- 



ial" periodicities along the rods and that its extension 
along the rods is fairly weak. For increasing Bjerrum 
length the above features get more pronounced and the 
oscillatory structure extends towards larger separations, 
while for decreasing Bjerrum length the features dimin- 
ish. At £b — ^ cr, at which the total force is approxi- 
mately 0, the only remaining feature of all three correla- 
tion functions is the correlation hole at 2 = (data not 
shown), which survives at even smaller values for £b- 



We would like to emphasize that this correlation hole 
at z = is not sufficient to entail attractive forces. This 
is important for the following reason: In the Wigner crys- 
tal approach the correlation energy of an ion can be well 
approximated by its interaction energy with the counter- 
charge within its correlation hole, which turns out to give 
the dominant contribution to the energy up to tempera- 

This may give 



tures far beyond the crystallization |52 



the impression that the existence of this first correla- 
tion hole is also sufficient for the appearance of attractive 
forces. Fig. |lj gives another illustration that this is not 
the case. The correlation functions gi and 172 are plot- 
ted for system 2 for three values of the Bjerrum length 
^b/c = 0.125, 0.5 and 3. As can be seen from Fig. ^, the 
first case corresponds to the strongest total repulsion, the 
second to an approximately vanishing net force and the 
third case to attraction. While the correlation functions 
in the first two cases are very similar, in the third case 
pronounced interlocking is observed. 



One word of caution should be addressed to the nu- 
merous theories which apply the interlocking picture to 
DNA-like systems. Our simulations show that the net 
attractive force in System 3 around ^b/c « 4 — 5 is only 
weak compared to the other two systems, and thus an 
interlocking pattern at theses values of £b is basically 
invisible. 
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V. SUMMARY, QUESTIONS AND OUTLOOK 

We have revisited the phenomenon of attractive in- 
teractions in systems of charged rods in the presence of 
muhivalent counterions. Within bulk systems this ef- 
fect is seen as a negative osmotic pressure. Separating 
between (i) the electrostatic and (ii) the entropic plus 
excluded volume component unveils the key role played 
by the short range repulsive forces. We further explored 
these components by measuring their contribution to the 
force between two parallel charged rods as a function of 
Bjerrum length £b- An important result is that an over- 
all attraction can be due to (i) both electrostatic and 
excluded volume forces bringing the rods together, (ii) 
an electrostatic attraction overcoming excluded volume 
repulsion, and (iii) a short range excluded volume force 
pushing the rods together against an electrostatic repul- 
sion. The excluded volume induced attraction can occur 
even if the average density between the rods is higher 
than outside, hence it is not a simple depletion mech- 
anism. This suggests that there exists more than one 
mechanism for attraction between rods. It is the subtle 
interplay between several effects which ultimately deter- 
mine the total force. We have shown that this crucially 
depends on the parameters of the system, like the line 
charge density or the rod radius. 

We also presented measurements of various ionic corre- 
lation functions. For the bulk systems the scaling of the 
first peak of the three-dimensional pair correlation func- 
tion indicates the existence of a three-dimensional corre- 
lated liquid, but only if one corrects for trivial effects due 
to the inhomogeneous average ion density. The assump- 
tion of a correlated liquid forming on the background of 
the rods is thus too simple. If the attractive interactions 
are still to be explained on the basis of the weak and 
masked three-dimensional correlations, this would imply 
that the Wigner crystal idea remains useful far beyond 
the ground state. 

For the two rod case we studied the distribution of ions 
around one rod relative to the other. We always found a 
correlation peak towards the other rod, but the rest of the 
distribution depends sensitively on Bjerrum length, line 
charge density and the rod radius. On the far side of the 
rod there exists a correlation hole at low Bjerrum length, 
which may turn into a peak at increasing Bjerrum length. 



Also, secondary peaks can appear at a certain angle with 
respect to the rod-rod direction. The periodicity of the 
azimuthal correlation function implies commensurability 
constraints, which influence the structure. The azimuthal 
correlation function was also shown to depend on the 
condensation radius 5r. It has been demonstrated for one 
system that the peak at ip = can turn into a correlation 
hole if one focuses on the ions most close to the rod. 
While this is important for the direction of the short 
range repulsive forces, it is presently unclear under which 
circumstances this effect occurs. 

We finally confirmed the recently proposed picture of 
ions mutually interlocking between the rods. The clos- 
est ions along the direction of the rod have a tendency 
to be condensed on different rods, even though defects 
in this structure occur quite frequently. In agreement 
with earlier observations on discretized models |56 this 
structure is remarkably short ranged. While attractive 
interactions have been found for all cases in which a first 
correlation peak is discernible, the mere existence of a 
first correlation hole was demonstrated to be insufficient. 

The phenomena observed in our molecular dynamics 
simulations also give rise to several new questions, for in- 
stance: How can one predict the relative contribution of 
electrostatic and short range forces to the net interaction 
between rods? The structural changes in the azimuthal 
correlation functions result from usual pair correlation 
functions being "wound up" around the rod. In conse- 
quence, the condensed ions predominantly push or pull 
from specific directions, not necessarily aligned with the 
direction joining the rods. What effects on the net force 
does this have? Under which circumstances do the ions 
between the rods get pulled away from the surface, and 
when is this the dominant force of attraction? Further 
studies along these lines are currently under way. 
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